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Distance  Calculation  Between  a Point 
and  a NURBS  Surface 


Eva  Dyllong  and  Wolfram  Luther 


Abstract.  In  this  paper,  we  consider  the  computation  of  an  Euclidean 
shortest  path  between  a point  and  a modelled  curve  or  surface  in  three- 
dimensional  space,  which  is  one  of  the  fundamental  problems  in  robotics 
and  many  other  areas.  A new  accurate  algorithm  for  the  distance-calcula- 
tion between  a point  and  a NURBS  curve  and  its  extension  to  the  case 
of  a point  and  a NURBS  surface  is  presented.  The  algorithm  consists  of 
two  steps,  and  is  crucially  based  on  appropriate  projections  and  subdivi- 
sion techniques.  To  solve  a nonlinear  polynomial  system  derived  from  the 
classical  formulation  of  the  distance  problem,  the  well-known  Newton-type 
algorithms  or  subdivision-based  techniques  first  considered  by  Sherbrooke 
and  Patrikalakis  are  used.  Their  modifications  in  conjunction  with  a low 
subdivision  depth  in  the  presented  algorithms  yield  a verified  enclosure  of 
the  solution. 


§1.  Introduction 

The  distance-calculation  is  an  essential  component  of  robot  motion  planning 
and  control  to  steer  the  robot  away  from  its  surrounding  obstacles  or  to  work 
on  a target  surface.  The  obstacles  may  be  polyhedral  objects,  quadratic  sur- 
faces, which  include  spherical  and  cylindrical  surfaces  or  more  general  surface 
types  like  the  non-uniform  rational  B-splines  (NURBS).  Most  of  the  well- 
known  algorithms  in  the  fields  of  computational  geometry,  robotics  and  Com- 
puter Aided  Design  are  focused  on  computing  the  distance  between  polyhedra, 
as  the  problem  is  easier  to  solve  and  the  answer  is  sufficient  for  many  prob- 
lems. For  example,  if  a free-form  designed  obstacle  like  a NURBS  surface  is 
located  at  a great  distance  from  a moving  robot,  then  it  is  sufficient  in  the 
next  step  to  know  the  distance  values  from  certain  sensor  points  on  the  robot 
to  the  convex  hull  of  the  NURBS  control  points,  which  describes  a convex 
polyhedron,  instead  of  the  more  time-consuming  and  expensive  computation 
of  the  exact  distance  values.  But  if  the  robot  approaches  an  obstacle,  more 
details  are  necessary,  and  fast  and  accurate  algorithms  for  finding  the  nearest 
point  on  the  NURBS  curve  or  surface  are  highly  recommended. 
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There  is  an  abundance  of  literature  to  calculate  the  distance  between 
convex  and  non-convex  objects.  For  convex  polyhedra,  a lot  of  algorithms 
exist  [1,4],  Two  well-known  algorithms,  the  Gilbert  method  (GJK,  [8])  and 
the  algorithm  of  Canny  and  Lin  (CL,  [3])  present  iterative  solutions  to  the 
problem,  i.e.,  both  construct  a sequence  of  pairs  of  proposed  distance  points 
which  are  then  improved  by  gradient  descent.  Another  wide  field  consists  of 
algorithms  for  general,  mainly  convex  objects  [5,7,13].  But  in  particular  for 
objects  defined  by  NURBS,  there  are  only  a few  contributions  in  the  litera- 
ture. In  [2],  Cameron  and  Turnbull  focus  on  computing  the  distance  between 
convex  objects  defined  by  NURBS  curves  or  patches,  for  which  the  critical 
step  is  the  evaluation  of  the  support  mapping.  The  method  is  based  on  the 
Gilbert-Foo  algorithm  [7]  for  general  convex  objects  with  adjustments  of  the 
termination  criteria  like  the  support  mapping  for  NURBS,  which  uses  deriva- 
tives of  their  basis  functions  and  the  Newton-Raphson  method  (NR  solver) 
for  finding  the  roots.  An  algorithm  for  the  computation  of  stationary  points 
of  a squared  distance  function  is  presented  in  [11].  This  problem  is  converted 
to  rie  polynomial  equations  with  variables  expressed  in  a tensor  product 
Bernstein  basis.  The  solution  method  uses  subdivision  relying  on  the  convex 
hull  property  of  Bernstein  polynomials  and  minimization  techniques. 

In  this  paper  we  describe  a new  algorithm  which  consists  of  two  steps,  and 
is  mainly  based  on  accurate  projections  and  subdivision  techniques.  In  the  first 
step,  the  NURBS  curve  is  decomposed  into  rational  Bezier  segments.  Then 
some  evaluations  of  suitable  scalar  products  decide  on  further  subdivision  of 
a rational  Bezier  segment.  This  subdivision  is  iterated  until  certain  criteria 
are  fulfilled.  In  addition,  a composition  of  the  new  method  together  with 
the  classical  formulation  of  the  distance  problem  based  on  the  calculation  of 
a solution  of  nonlinear  polynomial  systems  is  presented.  The  algorithm  is 
extended  to  the  case  of  a NURBS  surface. 


§2.  Problem  Formulation 

A NURBS  curve  C{u)  of  degree  p is  a vector-valued  function  of  one  parameter 
defined  by 


a < u < b. 


where  {Pi  = {xi,yi,Zi)  S R^}"_o  are  the  control  points,  {wjj-Lg  are  the 
weights,  and  {A^:,p(m)}"=o  are  the  pth-degree  B-spline  basis  functions  defined 
on  the  knot  sequence  U = with  {w;  = and  {uj  = 

Let  = (ujx,ujy,ijjz,uj)  = {x',y',z',uj')  € IR'*  and  H be  the  perspective  map 
given  by 

P ^ H{P-}  = H{(x\y',z',uj')}  = ( 

{{x',y',z'),  ifo;'=0. 

Applying  H to  the  nonrational  B-spline  curve  in  homogeneous  coordinates 

&{u)  = j2NiAu)Pr  (1) 

»=o 
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yields  the  corresponding  NURBS  curve  C{u),  i.e.,  ^(u)  = H{C‘^{u)}. 

Similarly  to  the  curve  case,  we  define  a NURBS  surface  using  the  tensor 
product  scheme.  Accordingly,  a NURBS  surface  S{u,v)  is  a bivariate  vector- 
valued piecewise  rational  function  of  the  form 

Si=o  !Cj=o  ^ 

E”=o  Er=o  NMNMwij  ’ ^ 

with  the  bidirectional  control  net  Pi,D  the  weights  Uij  and  the  B-spline  basis 
functions  iVj,p(u)  and  Nj^g{v). 

For  a given  point  ^ € IR®  we  jiddress  the  problem  of  finding  a shortest 
straight  line  segment  [Q,D\  with  D € C{u)  or  D E S{u,v),  respectively.  We 
assume  that  all  weights  of  the  rational  curves  and  surfaces  are  positive,  to 
ensure  that  the  convex  hull  property  holds. 

§3.  Distance  Algorithm  for  a Point  and  a NURBS  Curve 

In  this  section,  an  efficient  and  accurate  algorithm  for  distance  calculation 
between  a given  point  Q and  a NURBS  curve  C{u)  is  presented,  which  is  a 
kind  of  an  adaptive  system  of  solution  methods.  The  extension  to  the  case  of 
a NURBS  surface  works  analogously. 

The  algorithm  consists  of  two  steps.  In  the  first  step,  the  NURBS  curve 
is  decomposed  into  rational  Bezier  segments  Cj{u),  j = 1, . . . ,np,  which  can 
be  realized  once  in  the  preparation  phase.  In  [12]  Piegl  and  Tiller  present 
an  efficient  algorithm  for  computing  the  Up  Bezier  segments  using  the  ho- 
mogeneous form  given  by  (1).  Thus,  oiu-  task  is  to  compute  the  distance 
between  a point  and  a rational  Bezier  curve.  After  decomposition,  the  dis- 
tances between  Q and  each  endpoint  of  the  rational  Bezier  segments  Cj{u) 
are  calculated,  and  the  smallest  value  is  stored  as  a first  rough  approximation 
to  the  distance  value  d.  In  the  second  step,  the  rational  Bezier  segments  are 
processed  gradually.  Let  fc  = 0, . . . ,p,  be  the  control  points  of  the  j-th 
Bezier  segment  Cj{u),  j E {1, . . • ,rip}.  Then,  for  each  Pj,k,  fc  = 1, . . . ,p  — 1, 
the  distance  to  the  straight  line  supporting  the  line  segment  l{Pj,o,Pj,p)  be- 
tween the  endpoints  Pj^  and  Pj^p  is  calculated,  and  for  each  projection  point 
Rj,k,  fc  = 1, . . . ,p  — 1,  on  the  line,  we  test  if  Rj^k  belongs  to  the  line  segment 
l{Pjfi,Pj^p),  using  suitable  scalar  product  evaluations.  If  Rj^k  ^ KPjfijPj,p) 
for  at  least  one  fc  G {1, . . . ,p  — 1},  then  the  Bezier  segment  is  subdivided  into 
two  Bezier  segments,  for  which  the  second  step  of  the  algorithm  has  to  be 
started  again.  Otherwise,  the  algorithm  tests  whether  the  Bezier  segment  can 
be  replaced  by  the  line  segment  l{Pjfi^  Pj,p)  using  the  theorem  by  Wang  and 
Xu  (see  Sec.  3.2).  If  Cj{u)  is  nearly  a straight  line,  with  a given  accuracy 
e,  then  the  distance  between  the  point  Q and  the  line  segment  l{PjfijPj,p) 
is  calculated.  The  distance  d is  updated  if  a smaller  value  is  found,  and  j is 
replaced  by  j -|- 1.  If  Cj{u)  is  not  nearly  a line  segment,  the  following  scalar 
products 


S{u,v)  = 


Sk  ■=  {Pj,i  - Rj,i)  ■ {Pj,k  - Rj,i)  for  fc  = 2,...,p-l. 
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Fig.  1.  Convergence  of  the  algorithm. 


are  calculated.  If  all  Sk  > 0,  i.e.,  all  A:  = 1, . . . ,p  — 1,  are  on  the  same 
side  of  the  line  supporting  Pj,p),  then  the  algorithm  tests  the  position  of 

point  Q.  Otherwise,  the  Bezier  segment  Cj{u)  is  subdivided  into  two  segments. 
To  test  whether  Q lies  in  an  influence  area  of  the  Bezier  segment  Cj{u),  the 
projection  point  Rq  of  Q on  the  line  supporting  Pj,p)  is  calculated  and 

its  position  on  the  line  is  checked.  If  Rq  ^ l{Pjfi,Pj^p),  the  distances  d{Q,  Pj,o) 
and  d{Q,  Pj,p)  are  calculated,  d is  updated  if  necessary  and  j is  replaced  by 
j + 1.  Otherwise,  the  Bezier  segment  has  to  be  subdivided  into  two  segments 
to  increase  the  accuracy  of  the  result.  Fig.  1 shows  how  the  algorithm  works. 

The  subdivision  of  a rational  Bezier  segment  can  be  continued  until  the 
termination  criteria  (see  Sec.  3.2)  are  fulfilled  or  it  can  be  interrupted  after 
some  steps,  and  afterwards  the  distance  problem  can  be  transformed  in  terms 
of  the  solution  of  the  polynomial  equation 

(0  - C{u))  ■ C'{u)  = 0 (2) 

in  the  variable  u,  where  C'(«)  describes  the  derivative  of  the  curve  C[u).  This 
is  mainly  recommended  for  NURBS  surfaces  with  a large  curvature  to  avoid  a 
high  depth  of  subdivision.  We  calculate  the  roots  of  this  equation  using  either 
the  well-known  (interval)  Newton  method,  or  one  of  the  recently  implemented 
solution  methods  briefly  described  in  Sec.  3.3  (see  [10]).  A diagram  illustrating 
the  outline  of  the  distance  algorithms  is  given  in  Fig.  2. 


input:  ^(u):  n,p,  U, 


decomposition:  C(u)  l->  p,  Pjk,a)j j,  (j 


first  approach  to  distance  d :=  min(^  min  d(^,  „),  d(j^,  ^ ^ 

sub :-  0 ' 


for  7= 


yes 

~ ___jub  < __ 

no 

estimation  criteria 

subdivision  of  ?(u) 
sub++ 

NP  / LP  / PP  solver 

distance  calculation  to  line  segment 

update  d and  list  of  distance  points 

output:  distance  value  d,  list  of  distance  points 


Fig.  2.  Outline  of  the  distance  algorithm. 
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3.1.  Subdivision 

The  subdivision  of  the  control  polygon  of  C^(u)  determinated  by  {Pj^k}k=o 

into  two  Bezier  segments  with  control  points  {Qoi}f=o  {Qi  i}f=o  works 
in  homogeneous  coordinates,  and  reads  as  follows; 

for  k = 0 to  p do 

begin  Qo,k  ■=  PlM,p-k  :=  Pj,p-k\ 

for  1 = 0 to  p-k-1  do 

Pf.i  ■■=  (Pr,i  + ^;r/+i)/2'0  end 

Extending  the  idea  from  Bezier  curve  C^{u)  to  Bezier  surface  Sf  j{u,v)  by 
calling  the  routine  twice,  first  in  u and  then  in  v direction,  a subdivision  of 
Sf  j{u,v)  into  four  Bezier  patches  is  realized  (see  [12]). 

3.2.  Termination  criteria 

In  [15]  Wang  and  Xu  prove  the  following  theorem: 

Theorem  1.  For  the  rational  Bezier  curve  C(u)  of  degree  p, 

d{C{u),l{Po,Pp))  < $(wo,...,Wp)-  rnax  d{Pi,l{Po,Pp)), 

l<»<p—  1 

where  C(u)  = with  the  Bernstein  polyno- 

mials Bi^p(u)  of  degree  p,  and 

$(wo,...,Wp)  :=  1-  ^l  + max(wo\wp^)(^^max_^Wj)(2P~^  - 1)^  , 

d{C{u),l{Po,Pp))  :=  sup  { inf  d((?(^t),^Po  + (l-<)-Pp)}■ 

a<u<6 

If  after  some  subdivision  steps,  d{Q,l{Pjfi,Pj^p))  >d  + d{Cj{u),l{Pjfi,Pj^p)), 
or  the  curve  can  be  approximated  by  l{Pjfi,Pj^p),  i.e.,  d{Cj{u),l{Pj^o,  Pj^p))  is 
not  greater  than  the  desired  tolerance  e,  the  subdivision  of  the  segment  stops 
after  this  step. 

If  LOi  = const,  the  curve  Cj{u)  describes  a Bezier  curve,  and  the  termi- 
nation criterion  of  Theorem  1 is  reduced  to  testing  the  following  conditions: 
(1  — l/2P~^)d{Pj^kP{Pj,o,Pj,p))  < e for  fc  = 1, . . . ,p  — 1.  In  addition,  for  a 
Bezier  curve  the  following  theorem  proved  in  [14]  specifies  the  number  of  nec- 
essary subdivisions,  i.e.,  after  r subdivision  steps  the  curve  can  be  replaced 
by  the  line  segments: 

Theorem  2.  For  the  Bezier  curve  C{u)  of  degree  p with  control  points 
{Pi  = (xi,yi,Zi)}f^g  and  any  given  e > 0,  let  L :=  maxo<i<p_2{|a:i  - 2xi+i  -1- 
Xi+2\,  \yi-2yi+i+yi+2\,  \zi-2zi+i+Zi+2\},  and  r :=  \og^{y^p{p-l)L/{8e)). 
Then,  for  a < a < P < b and  — log2((/8  — «)/(6  — a))  > r, 

d(C(u),liC{a),C(m<£- 

If  at  the  beginning  of  the  subdivision  the  value  r can  be  calculated,  then  it  is 
not  necessary  to  perform  the  termination  tests  of  Theorem  1. 
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3.3  Solution  techniques  and  complexity  analysis 

The  distance  problem  can  be  converted  into  a problem  of  computing  all  roots 
of  a system  of  nonlinear  polynomial  equations  in  one  or  two  variables.  There 
are  two  techniques  designed  to  solve  such  a problem  in  variables  efficiently; 
the  projected-polyhedron  (PP)  and  the  linear  programming  (LP)  technique, 
developed  by  Sherbrooke  and  Patrikalakis  [10].  They  rely  on  representation 
of  polynomials  in  the  multivariate  Bernstein  basis,  the  convex  hull  property 
and  on  the  subdivision  or  linear  programming  technique.  Alternatively,  the 
Newton-Raphson  method  can  be  used  to  find  the  roots  of  the  nonlinear  poly- 
nomial equations. 

Next,  we  analyse  the  amount  of  time  required  to  execute  each  step  of 
the  distance  algorithm  in  Ccise  of  the  NURBS  curve.  The  decomposition  of 
the  curve  C{u)  of  degree  p into  Up  Bezier  segments  takes  at  most  0{p  ■ Up) 
operations,  and  the  first  approach  to  the  distance  value  needs  Up  + \ steps.  The 
total  cost  of  the  distance  calculation  for  rip  Bezier  segments  with  subdivision 
depth  of  k is  in  worst  case  0(2*p^  • rip)  independently  of  the  method  used 
(subdivision-based  technique  or  PP/LP  solver)  for  finding  the  distance  points. 
In  this  case  (rig  = 1),  close  to  a simple  root,  quadratic  convergence  is  achieved. 


§4.  Distance  Algorithms  for  a Point  and  a NURBS  Surface 

In  the  case  of  a NURBS  surface  S{u,v),  the  distance  algorithm  maintains  its 
structure.  After  the  decomposition  of  the  surface  into  rip  • Uq  Bezier  patches 
Sij{u,v)  with  control  points  {0  < k < p,  0 < I < q),  the  first  approxi- 
mation to  the  distance  value  d is  calculated  taking  the  minimum  of  distances 
d{Q,  Pk\)  for  k = 0,p  and  I = 0,q.  The  termination  critera  for  subdivision  of 

a non-degenerate  Bezier  patch  Sij{u,v)  {Pqq,  Pp{,,  pQq  are  not  collinear)  are 
modified  in  the  following  way: 


and  {P^  - P^i) . ((p;;j  - P';j)  X (P';j  - P’;j)) 


< e 


for  all  1 < A:  < p,  1 < / < g,  and  e = 0,p  (where  x denotes  the  cross 
product).  In  this  case  Sij{u,v)  can  be  replaced  by  a plane  segment  defined 
by  Pg  0,  Pp’o,  and  Pg’^,  the  subdivision  of  Sij{u,  v)  is  stopped,  and  the  distance 
between  the  plane  segment  and  the  point  Q is  calculated.  The  remaining  tests 
are  performed  in  the  u and  v directions  analogously  to  the  curve  case  using 
the  tensor-product  structure  of  Si,j{u,v).  If  the  point  Q does  not  lie  in  the 
influence  area  of  the  Bezier  patch  Sij{u,v),  i.e.,  the  projections  onto  the  lines 
forming  the  boundary  of  the  triangle  defined  by  Po'g,  Pp  l,  and  Pg’^  or  defined 
by  Pp’jj,  Pp  Q and  Pg’^  do  not  belong  to  sides  of  the  triangle,  the  distance 
between  the  point  Q and  one  of  the  boundary  lines  Sij{u,v),  u € {a,  6}, 
V e {c,  d},  is  calculated,  d is  updated  if  necessary,  and  the  subdivision  of 
Sij{u,v)  is  interrupted. 
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If  a particular  depth  of  subdivision  is  obtained,  the  subdivision  of  the 
Bezier  patch  can  be  stopped,  and  the  PP  or  LP  solver  for  nonlinear  polynomial 
systems  in  two  variables  can  be  applied  to  the  Bezier  patch.  The  equations 
for  Sij{u,v)  read  as  follows  analogously  to  (2): 


2p-l  2q  2p  2q-l 

'^0'k,iBk,2p-i{u)Bi2q{v)  = 0 and  ^ ^ h,iBkfip{u)Bi^2q-i{v)  = 0 


fc=0  1=0 


fc=0  J=0 


with 

min(p— l,fc) 

«fc.i  = XI 

5=max(0,fc— p)t=max(0,/— g) 
min(p,fe)  m\n{q—l,l) 

E E - <?)■ 

s=max(0,k—p)t=mSLx(0,l—q)  \k)\  I ) 


min(9,()  ('P-1')  f f ® '1 

E\  3 ) \k-a) / fji,j  _ pi,j\/pi,j 


The  combination  of  a classical  formulation  of  the  distance  problem  and 
the  subdivision  technique  is  recommended  if  a high  subdivision  depth  is  ex- 
pected, e.g.,  if  r in  Theorem  2 is  too  large  in  case  of  a very  bent  Bezier  curve. 


§5.  Concluding  Remarks 

The  method  described  in  this  paper  computes  the  distance  between  a point 
and  a NURBS  curve  or  surface.  Our  goal  was  to  provide  a reliable  method  to 
solve  this  problem.  The  first  few  subdivision  steps  and  tests  quickly  locate  the 
regions  of  potential  solutions.  Then  the  subdivision  can  be  either  continued, 
or  one  of  the  equation  solvers  or  even  a distance-calculation  algorithm  for 
polyhedra  can  be  applied  [4,5].  We  have  developed  an  interval  version  of  the 
PP/LP  algorithm  using  interval  arithmetic  and  considering  a correct  handling 
of  roots  of  order  two,  suitable  modifications  of  Graham’s  scan  algorithm  for 
building  the  convex  hull,  the  revised  simplex  method  by  Gass,  and  an  adapted 
interval-based  subdivision  by  de  Casteljau.  The  solver  has  been  implemented 
in  C-I--I-  using  the  library  Profil/BIAS  (see  [9]).  This  improves  the  robustness 
of  the  distance-algorithm,  assures  an  interval  enclosure  of  the  solution,  and 
makes  it  suitable  for  verification  of  off-line  tasks  in  path  planning. 

Some  modifications  to  the  algorithms  could  improve  performance,  e.g.,  if 
upper  bounds  on  the  derivatives  of  order  two  for  the  curve  or  surface  are  known 
[6].  But  doubtless  our  NURBS-based  algorithm  will  be  slower,  e.g.,  compared 
with  our  algorithms  [4,5],  where  we  deal  with  the  objects  as  polyhedra.  In  a 
more  complete  distance  tracking  system,  such  as  a manipulator  in  a complex 
environment,  a progressive  switch  from  a spherical  or  polyhedral  enclosure  of 
the  objects  to  NURBS  surfaces  is  recommended,  especially  if  contact  problems 
are  investigated. 
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